function [P, y, x, MVndx] = VarTauchen(Ns,F0,F,bandwidth,Sigma)

tol = 10^-5;

% y_t = F0 + F*y_{t-1} + \varepsilon

%% construct the range over which each of the compents may vary
% Y_t = F0 + F Y_{t-1} + Q \varepsilon_t
Ny       = size(F, 1);
Nstar    = Ns^Ny;
Q        = chol(Sigma)'; % QQ' = Sigma
iQ       = inv(Q);
Iy       = eye(Ny);
EY       = (Iy - F)\F0;
VarY = uncvary(F,Sigma);

% X_t = Q^{-1} Y_t = F0x + Fx X_{t-1} + \varepsilon_t 
Fx       = iQ * F * Q;
F0x      = iQ * F0;
EX       = iQ * EY;
VarX     = iQ * VarY * iQ'; % = dlyap(F, Iy);
StdX     = sqrt(diag(VarX));



%% Define Grids:
Ny = size(F,1);

Ub = EX + bandwidth.*StdX;
Lb = EX - bandwidth.*StdX;

steps   = (Ub - Lb) / (Ns - 1);

griddy  = NaN(Ns,Ny);
for x = 1 : Ny
   griddy(:,x) = Lb(x) : steps(x) : Ub(x);
end

%% Index for Multivariate Grid
MVndx                = MultIdx(Ns, Ny);

% note: midpoints also used for conditional means! (see below)
x = griddy(MVndx);
y = x * Q';

endpoints = griddy(1:end-1,:) + repmat(steps' / 2, Ns-1,1);


condmean = x * Fx' + repmat(F0x',Nstar,1);
condstd  = ones(1, Ny);
P = NaN(Nstar, Nstar);
for s = 1 : Nstar
   probby   = diff([zeros(1,Ny); ...
      normcdf(endpoints, repmat(condmean(s,:), Ns-1, 1), repmat(condstd,Ns-1,1)); ...
      ones(1,Ny)]);
   P(s,:)   = prod(probby(MVndx), 2)';
end



%GridY = 


    function MIdx = MultIdx(Ns,Ny)
        Nstar = Ns^Ny;
        MIdx = NaN(Nstar,Ny);
        UVndx = (1:Ns)';
        MIdx(:,1) = repmat(UVndx, Ns^(Ny-1), 1);

        for j = 2 : Ny
           n1          = Ns ^ (j - 1);
           n2          = Ns ^ (Ny - j);
           MIdx(:,j)  = repmat(kron(UVndx, ones(n1,1)), n2, 1) + (j-1) * Ns;
        end
    end

    function VarY = uncvary(F,Sigma)
    % This subfunction computes the unconditional variance of Y using the
    % iteration suggested in Tauchen:

    % VarY(r) = F VarY(r-1) F' + Sigma
    err = 1; OldVarY = eye(size(Sigma));
    while err > tol
        NewVarY = F*OldVarY*F' + Sigma;
        err = max(max(abs(OldVarY-NewVarY)));
        OldVarY = NewVarY;
    end
    VarY = NewVarY;
    end

end
    



